rm(list=ls())
library('GJRM')
load("dataBJPOLS.RData")

last.cases$vote_post <- as.numeric(last.cases$vote_post > 0.80)
last.cases$vote_pre <- as.numeric(last.cases$vote_pre > 0.80)
last.cases$regis_before <- as.numeric(last.cases$regis_before > 0.80)

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"
outc.1 <- "vote_post"

out.eq <- vote_post ~ pti +  as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)  + age + I(age^2) +  
  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before +
  as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)                      

treat.eq <- pti ~ judgeiv_hd +  as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity) + age + I(age^2) +  
  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before +
  as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)

f.list <- list(treat.eq, out.eq)

mr <- c("probit", "probit")
bvp <- gjrm(f.list, data = last.cases, model ="B", margins = mr)

last.cases_1 <- last.cases_2 <- last.cases
last.cases_1$pti <- 1
last.cases_2$pti <- 0

out1 <- predict(bvp, eq = 2, newdata = last.cases_1, type = 'response')
out0 <- predict(bvp, eq = 2, newdata = last.cases_2, type = 'response')
outF <- mean(out1) - mean(out0)
outF

sims <- 500
library(parallel)
library(doParallel) 
cl <- detectCores()/2

registerDoParallel(cl) 

foreach(j = 1:sims) %dopar% {
  out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]

  out.eq <- vote_post ~ pti +  as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity) + age + I(age^2) +  
    as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before +
    as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)                      
  
  treat.eq <- pti ~ judgeiv_hd +  as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity) + age + I(age^2) +  
    as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before +
    as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)
  
  f.list <- list(treat.eq, out.eq)
  
  mr <- c("probit", "probit")
  bvp <- gjrm(f.list, data = last.cases.boot, model ="B", margins = mr)
  
  last.cases_1 <- last.cases_2 <- last.cases.boot
  last.cases_1$pti <- 1
  last.cases_2$pti <- 0
  
  out1 <- predict(bvp, eq = 2, newdata = last.cases_1, type = 'response')
  out0 <- predict(bvp, eq = 2, newdata = last.cases_2, type = 'response')
  out <- mean(out1) - mean(out0)
  save(out, file = paste0("./BiProbitResults/BiProbitRes_", round(runif(1),10) * 1e10, ".RData"))
}

stopImplicitCluster()

library('GJRM')
load("dataBJPOLS.RData")

last.cases$vote_post <- as.numeric(last.cases$vote_post > 0.80)
last.cases$vote_pre <- as.numeric(last.cases$vote_pre > 0.80)
last.cases$regis_before <- as.numeric(last.cases$regis_before > 0.80)

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"
outc.1 <- "vote_post"

out.eq <- vote_post ~ pti + as.factor(court_time1_ds) + as.factor(court_time2_ds) + age + I(age^2) +  
  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before +
  as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)                      

treat.eq <- pti ~ judgeiv_hd +  as.factor(court_time1_ds) + as.factor(court_time2_ds)  + age + I(age^2) +  
  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before +
  as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)

f.list <- list(treat.eq, out.eq)

mr <- c("probit", "probit")
bvp <- gjrm(f.list, data = last.cases, model ="B", margins = mr)

last.cases_1 <- last.cases_2 <- last.cases
last.cases_1$pti <- 1
last.cases_2$pti <- 0

out1 <- predict(bvp, eq = 2, newdata = last.cases_1, type = 'response')
out0 <- predict(bvp, eq = 2, newdata = last.cases_2, type = 'response')
outF <- mean(out1) - mean(out0)
outF

list.of.files <- list.files("./BiProbitResults/")
results <- list()

for(i in 1:length(list.of.files)) {
  load(paste0("./BiProbitResults/", list.of.files[i]))
  results[[i]] <- out
}

se.bip <- sd(do.call('rbind', results))

outF
se.bip

